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Abstract. 

The first consistent model for the turbulent tachocline is presented, with the tur- 
bulent diffusivity computed within the model instead of being specified arbitrarily. 
For the origin of the 3D turbulence a new mechanism is proposed. Owing to the 
strongly stable stratification, the mean radial shear is stable, while the horizontal 
. shear is expected to drive predominantly horizontal, quasi-2D motions in thin slabs. 

Here I suggest that a major source of 3D overturning turbulent motions in the 
£f>l , tachocline is the secondary shear instability due to the strong, random vertical shear 

■ arising between the uncorrelated horizontal flows in neighbouring slabs. A formula 
■^j- ' for the vertical diffusivity due to this turbulence, equation (9), is derived and applied 
£f~) , in a simplified ID model of the tachocline. It is found that Maxwell stresses due to 

an oscillatory poloidal magnetic field of a few hundred gauss are able to confine 
the tachocline to a thickness below 5 Mm. The integral scale of the 3D overturning 
turbulence is the buoyancy scale, on the order of 10 km and its velocity amplitude 

■ is a few m/s, yielding a vertical turbulent diffusivity on the order of 10 8 cm 2 /s. 

. Keywords: Sun: interior, MHD, tachocline 
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I 1. Introduction 



In most current dynamo models (see e.g. Petrovay, 2000, for a re- 
view), the strong toroidal magnetic field responsible for solar activity 
$J< ■ is generated in the thin rotational shear layer below the convective 



zone (CZ), known as the tachocline. To hold sufficient magnetic flux 
to explain the flux observed at the surface without unrealistically high 
field strengths, this layer should be at least a few megameters thick. On 
the other hand, a magnetic field oscillating with a circular frequency 
io cyc = 2n/P, {P = 22 years is the solar cycle period) is known to 
penetrate a conductive medium only down to a skin depth of 

H skin = (2r ] /u; cyc ) 1 ^ (1) 

where rj is the magnetic diffusivity. If i/ s kin > 1 Mm then we conclude 
from this that for present-day dynamo models to work, we must have 
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7] ^ 10 8 cm 2 /s in the tachocline layer. For comparison, the molecular 
magnetic diffusivity is ~ 10 3 cm 2 /s in this part of the Sun, while the 
turbulent diffusivity in the convective zone is estimated to be 10 12 - 
10 13 cm 2 /s. 

Helioseismic evidence indicates that at low heliographic latitudes 
the middle of the tachocline is situated at 0.69 Rq (1 Rq = 698 Mm 
is the solar radius) and its thickness is less than a few percent of Rq 
(Kosovichev, 1996). There are indications that at high latitudes the 
tachocline is situated at slightly shallower depths, r = 0.705i? Q , and it 
may also be marginally thicker (Basu and Antia, 2001). Comparing this 
with the radius of the bottom of the adiabatically stratified convective 
zone, 0.71 Rq, it is clear that most of the tachocline lies in the radiative 
interior, especially near the equator, so the relatively high value of the 
(presumably turbulent) magnetic diffusivity inferred above is puzzling. 
Hence, there is clearly a need to understand the nature of small-scale 
turbulent motions in the tachocline. Unfortunately, while the thermal 
stratification in the tachocline is relatively well known, its fluid dy- 
namical properties, including the precise profile of the rotational flow 
v(r,6), the meridional flow, and the characteristics of turbulence, are 
very poorly constrained by observations. 

The tachocline is thus basically an MHD shear flow in the az- 
imuthal direction in a thin, rotating spherical shell with both radial 
and latitudinal shear and a strongly stable stratification (Richardson 
number ~ 10 3 ). Theoretical understanding and experimental knowl- 
edge regarding turbulence in stratified shear flows has been reviewed 
by Hopfinger (1987) and Thorpe (1987). These reviews ephasize the 
tendency towards two-dimensional, predominantly horizontal motions 
in such situations, the vertical motions being limited to the small 
buoyancy scale (see below). In line with the concepts outlined in those 
reviews, the present paper to relies on theoretical arguments concerning 
the properties of turbulent motions in such conditions. On the basis 
of such arguments, simplified one- or two-dimensional models may be 
constructed for the mean flow, or appropriate subgrid closure schemes 
may be constructed for full 3D large-eddy simulations. While such argu- 
ments have been attempted in the past (Spiegel and Zahn, 1992, Gough 
and Mclntyre, 1998), in the lack of detailed turbulence calculations to 
be used for subgrid closures, all tachocline models published to date 
have either simply ignored turbulence (Riidiger and Kitchatinov, 1997, 
MacGregor and Charbonneau, 1999, Garaud, 2001a, Cally, 2001), or 
assumed arbitrary fixed scalar turbulent diffusivities in 2D mean flow 
models (Forgacs-Dajka and Petrovay, 2001, 2002, Forgacs-Dajka, 2003) 
and in 3D LES (Miesch, 2001, 2002). 
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The aim of the present work is to attempt to remedy this situation 
by considering, on the basis of the known stability criteria, the possible 
sources of turbulence in a strongly stably stratified shear flow with both 
vertical and horizontal shear, and discussing the expected properties of 
the turbulence generated by it, on the basis of a dimensional analysis 
of the K—e equations. As our analysis does not consider the effects of 
spherical geometry, rotation, meridional circulation, non-adiabatic con- 
vective overshoot, or magnetic instabilities, it should only be regarded 
as a fist step towards a more comprehensive theoretical analysis of 
the problem of turbulence in the solar tachocline. These theoretical 
arguments can be found in Section 2. Then, in Section 3, as an illustra- 
tion of the use of such theoretical considerations, our prescription for 
calculating the turbulent diffusivity in the tachocline, equation (9), is 
incorporated in a simplified one-dimensional model for the tachocline. 
The results show that Maxwell stresses due to an oscillatory poloidal 
magnetic field of a few hundred gauss (a rather moderate value) are able 
to confine the tachocline to a thickness below 5 megameters. The con- 
ditions of validity of our treatment are discussed in Section 4. Finally, 
Section 5 concludes the paper by comparing our suggested picture of 
tachocline turbulence with numerical simulation results and discussing 
the implications of the results. 

2. A possible source of turbulence: the secondary shear 

instability 

The thermal stratification of the Sun is quite accurately known from 
a comparison of standard solar models with helioseismic inversion re- 
sults. The radiative interior below the CZ is characterized by significant 
negative values of AV. For rough estimates, a useful approximation in 
the tachocline region (i.e. near the top of the radiative zone) is AV ~ 
—0.015 2 [Mm], where z = r^ cz — r is the depth below the bottom of the 
convective zone, at a radius value of r^ cz = 0.71R Q . On the other hand, 
the pressure scale height in this region is Hp = —Pdz/dP ~ 50 Mm. 
With g = 5 • 10 4 cm 2 /s, this yields a Brunt- Vaisiila frequency 

iVgy [s~ 2 ] = -AVg/Hp ~ 1.5 • 10~ 7 z [Mm], (2) 

i.e. at a depth of a few megameters Abv ~ 10 _3 s _1 . 

A displaced mass element will then oscillate around its equilibrium 
position under the action of buoyancy on a timescale Agy ~ 10 3 s. The 
amplitude of the oscillation is clearly ~ v z /Nby, so in the presence 
of turbulent motions, these motions will be limited to a vertical scale 
lb = K 1 / 2 /N-qv, called the buoyancy scale. (K = v\ is the kinetic 
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energy in the vertical component of motions.) On the other hand, 
as an elementary estimate gives v ~ K 2 /e for the vertical turbulent 
diffusivity, one has If, ~ (y /N-qv) 1 / 2 . Vertical overturning motions on 
scales exceeding If, will be strongly damped by gravity wave emission. 

The pole-equator difference in the rotational rate of the convective 
zone is about 30 % of the equatorial rotation rate. Taking half of this 
value to be the characteristic amplitude of the differential rotation (cf. 
eq. (5) below), we have a differential rotation amplitude Avq ~ 3 • 
10 4 cm/s at the top of the tachocline. This value is clearly also roughly 
the amplitude of the vertical velocity difference across the tachocline, 
so the characteristic vertical shear is 

S ~ Av /H ~ 3 • 1(T 5 s- 1 (3) 

where H is the scale height of Av, as given by equation (5) below. 
Helioseismic calibrations (Basu and Antia, 2001) indicate H ~ 10 Mm. 
This yields a Richardson number Ri= N^y/S 2 ~ 10 3 . This enormous 
value shows that the vertical shear cannot directly drive turbulence in 
the tachocline. 

Buoyancy, however, cannot stabilize the horizontal shear. While 
linear stability analysis (Dziembowski and Kosovichev, 1987, Char- 
bonneau, Dikpati, and Gilman, 1999) shows that the horizontal shear 
is close to marginal stability in the 2D nonmagnetic case, nonlinear, 
3D effects and magnetic fields are known to lead to strong instability 
(Gilman and Dikpati, 2000, Dikpati and Gilman, 2001). The motions 
driven by the horizontal shear instability are predominantly horizontal, 
and their spatial scale is 

lh~Ro, (4) 

On the basis of general experience with shear flow turbulence 1 , the ve- 
locity scale Vh of the horizontal motions driven by the shear instability 
may be assumed to be close to the amplitude Av of the horizontal shear 
at the given depth: 

v h ~ Av = [uj(z, 6 = 0)- uj(z, 9 = ir/2)]R Q /2 (5) 

Similarly, their correlation time is assumed to be 

l h /v h ~ R Q /Av (6) 

Overturning turbulent motions in the vertical direction are impeded 
by the stable stratification, their scale being limited to If,. Owing to the 

1 In fact it is not obvious that this experience is a useful guide under conditions 
in the solar tachocline. This point will be discussed further in Sect. 4 below. 
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low vertical diffusivity and the finite correlation time, the horizontal 
motions will be characterized by a limited vertical correlation length 

lc ~ {vlh/vh) 1 ' 2 . (7) 

The random horizontal flows driven by the shear will then be limited 
to thin sheets of thickness l c , the motion in neighbouring sheets be- 
ing independent. This will give rise to random vertical shear between 
neighbouring sheets, of amplitude 

S 2 ~ vh/lc ~ {vl/vh) 1 / 2 (8) 

This secondary vertical shear is much stronger than the primary (mean) 
vertical shear, the corresponding Richardson number being 

Ri 2 = gAVul h /Hpvl 

Substituting here the characteristic values of the parameters, we find 
that Ri2 < 0.25, i.e. the secondary shear is unstable, if 

V< " a 4l h N* v - W z[Mm] W 

In a depth of a few megameters this value is u CI ~ 10 8 cm 2 /s, much 
higher than the molecular value, so we expect that, in the absence of 
other sources of turbulence, the secondary shear is strongly unstable. 

What is the characteristic amplitude of the turbulent motions driven 
by the secondary shear instability? In principle, this could be derived 
from a K—e model (or, more, generally, from a Reynolds stress model 
(Speziale, 1991). Assuming plane parallel geometry for simplicity, the 
relevant equations are of the general form 

Here, the non-local fluxes or third order moments are 

F K =^I (12) 

F e = vlel (13) 

ei being the local dissipation rate, while e is the mean dissipation. The 
production terms are usually modelled as 

Pk = -S\ (14) 
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P e = Cp-^Pk (15) 

while, assuming an ideal gas and the Boussinesq approximation p'/p = 
—T'/T (prime denotes fluctuations), the dissipation/sink terms read 

D K = e + g^T (16) 

D e = C D1 ^ + C D2 ^g^T (17) 

(The first terms on the r.h.s. represent viscous dissipation, while the 
second terms correspond to gravity wave emission. Note that in a suba- 
diabatic environment v z T' > 0, i.e. downmoving fluid parcels are hotter 
than average.) 

Performing a dimensional analysis on these equations we find that 
the diffusive timescale, corresponding to the non-local terms, is d 2 jv J> 
d 2 /v C r ~ 10 10 s, while the timescale associated with the shear produc- 
tion term is KjvS\ ~ Nbv/S% ~ N B vlhv/vl ~ lCT 5 ^ < 10 3 s. This 
implies that the transport terms can be neglected in the K-e equations. 
Under such conditions the equations have no stationary solution, as the 
values of the constants Cp,Coi and C02 ar e different in general. The 
intensity of turbulence will then keep increasing until the turbulent 
diffusivity reaches the critical value u CT , when further shear production 
is switched off. 

We thus conclude that the secondary vertical shear instability can 
be expected to drive overturning turbulence to the level v = v cr on a 
short timescale. The turbulence generated by this mechanism may then 
be crudely represented by the vertical diffusivity value given by eq. (9). 



3. Tachocline model 

3.1. Model equations 

We now proceed to develop a simple one-dimensional model for the so- 
lar tachocline, assuming that the secondary shear instability discussed 
in the previous section is the only source of turbulence in the tachocline 
region. Our computational domain will be restricted to the top of the 
radiative interior, below r\, cz . 

The convective zone is characterized by extremely high turbulent 
diffusivities. Due to the Coriolis force, the Reynolds stress tensor also 
has significant nondiagonal components in the convective envelope. 
These components imply an angular momentum transport which is 
thought to be the main driver of solar differential rotation. Based on 
the discussion of the previous section we expect that the amplitude 
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of turbulence in the tachocline is much lower than in the CZ. Thus, 
from the point of view of tachocline modelling, it is not unrealistic 
to regard the latitudinal differential rotation at r = r^ cz as a given 
boundary condition imposed at the top of the region of interest. This 
is tantamount to assuming that differential rotation is driven by a 
highly effective mechanism in the convective zone which is not seriously 
influenced by the processes in the tachocline. 

As the layer studied is thin, we also adopt a plane parallel rep- 
resentation for it, with constant density. (The effects of the stable 
density stratification are only implicitly taken into account by its role 
in determining the turbulent viscosity, eq. (9).) 

Thus, we regard the following model problem. Consider a plane 
parallel layer of incompressible fluid of density p, where the viscosity v 
and the magnetic diffusivity 77 depend on z only. At z = where z is 
the vertical coordinate (corresponding to depth in the solar application 
we have in mind) a periodic horizontal shearing flow is imposed in the 
y direction: 



(so that x will correspond to heliographic latitude, while y to the 
longitude). We assume a two-dimensional flow pattern (d/dy = 0) 
and v x = v z = (no "meridional flow"). An oscillatory horizontal 
"poloidal" field is is present in the x direction, given by 



(in Alfvenic units). The "toroidal" (i.e. y) component A of the vector 
potential obeys the corresponding component of the integral of the 
induction equation, which in our case simplifies to a diffusion equation: 



The upper boundary condition is A = Aocos(ut) at z = 0, where Aq 
fixes the poloidal field amplitude. 

The evolution of the azimuthal components of the velocity and the 
magnetic field is described by the corresponding components of the 
equations of motion and induction, respectively. Introducing v = v y 
and using Alfven speed units also for the toroidal magnetic field 



v y o = vq cos(kx) 



1 dA 
(^p) l / 2 ~dz~ 



(18) 




(19) 



b = B y {Airp)- l /\ 



(20) 



these can be written as 




(21) 
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db Jr ( ^dv d ( db\ 
01 

where we have taken into account that, owing to the thinness of the 
tachocline, the vertical derivatives dominate the diffusive terms. 2 As 
the imposed poloidal field V p is independent of x, Fourier transforming 
(21) and (22) in terms of x yields the same equations for the Fourier 
amplitudes v and b, except that d/dx is substituted by ik: 



- = ^y p cosM) + -^-j (23) 



db , Tr , . d ( db\ 

- = ^cosM) + -^-j (24) 

(Hats are omitted to simplify notation.) For a rough estimate, we write 
nb IP for the l.h.s. of (24), then substitute the resulting expression of 
b into (23), take the real part, and omit the factor cos 2 (out) in the first 
term on the r.h.s.: 

dv , 2t/2 „ . d ( dv\ 



» = -^ F » + M'&J (25) 

The equations to solve are thus (18), (19) and (25), with the turbu- 
lent diffusivities v = r] = v m + u CT given by equation (9) with the 
identification Vh = v. u m is a minimal diffusivity value ("molecular 
diffusivity" ) . 

The simplified form of the first term in equation (25) will not allow 
a correct reproduction of the periodic part of the time dependence. 
The important point is, however, that this sink term, representing 
the reduction of horizontal shear by Maxwell stresses, has the right 
amplitude and the correct scaling with V p , P, v, and k, so it may 
be expected to reproduce well the cycle- averaged flow amplitude as 
a function of z, which is our main interest here. Indeed, solving our 
equations for the case of constant diffusivities v and rj, the results are 
in a remarkably good agreement with the full 2D solutions in spherical 
geometry, presented in Forgacs-Dajka and Petrovay (2002). 

The equations were solved numerically by a finite difference scheme 
second-order accurate in time. All quantities were set to zero at the 
lower boundary, situated at zq = 60 or 30 Mm below r^, cz , while the 
boundary conditions applied at top (z = 0) were v = vo = 3 ■ 10 4 cm/s 
and different prescribed values for Aq. As equation (9) is singular at 
z = 0, v was set to a high finite value f max here. For v m we used 

The large-scale quasi- 2D motions are assumed not to contribute to angular 
momentum transport here, as the amplitude and sense of such a contribution is still 
highly controversial, cf. Spiegel and Zahn (1992), Gough and Mclntyre (1998). 
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2 [Mm] z [Mm] 

Figure 1. Horizontal differential rotation amplitude v — Vh, as denned in eq. (5) 
(top row), vertical turbulent diffusivity v (middle row), and poloidal magnetic field 
in Alfvenic units (bottom row), averaged over a solar half-cycle, as functions of 
depth below the convective zone, for a very low (left-hand column) and a medium 
(right-hand column) value of the field strength imposed at the top. Note that by 
coincidence, in the solar tachocline the field strength in gauss is roughly equal to 
the corresponding Alfven speed in cm/s. 
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5 10 15 20 25 30 5 10 15 20 25 30 

z [Mm] z [Mm] 



Figure 2. Same as Fig. 1 for two higher field strength values. 



the value 10 7 cm 2 /s. This is much higher than the actual molecular 
diffusivities in the tachocline, but using a realistic value would lead to 
forbiddingly long integration times. Similarly, a too high value for f max 
would lead to very short timesteps, also increasing the computing time 
to unaffordable values. Test runs with varying values of z^max 

and u m , 

however, show that these choices do not significantly distort the results. 
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Starting from an arbitrary intial state, the system was allowed to 
evolve until very nearly periodic behaviour sets in (in about 10 years, 
depending on the value of f m ), then average quatities for one 11-year 
half-cycle were computed and plotted as functions of depth (Figs. 1 
and 2). 

3.2. Results 

In the case with a very weak magnetic field (left-hand column in Fig. 1), 
it is straightforward to show that equation (25) with v given by (9) 
admits the analytic solution 



confirmed by the numerical calculations. This essentially means that 
in this case the shear penetrates as far down into the radiative interior 
as the placement of the lower boundary condition allows. The weak 
magnetic field itself only penetrates down to the skin depth given by 
equation (1), as expected. It is worth noting that the poloidal field 
amplitude shows a non-monotonic behaviour with depth in all cases, 
reaching its maximum at some finite z value. This is due to the vari- 
able diffusivity: the horizontal field lines tend to "pile up" where the 
diffusivity is significantly reduced. 

From the right-hand column of Fig. 2 we can see that a poloidal 
magnetic field of a few hundred gauss (peak strength 500 G) can confine 
the tachocline to a thickness of barely 5 Mm. 

One might think that an intermediate field strength might lead to a 
somewhat thicker tachocline. This is, however, not the case: an inspec- 
tion of the full series of results in our figures clearly shows that a weaker 
field simply results in an "aborted tachocline" , i.e. the horizontal shear 
is first reduced by a factor depending on the field strength in a thin 
layer of a few Mm, but below that layer, as the magnetic field is damped 
by the skin effect, it follows the field- free solution (26), with a lower 
amplitude. 




(26) 



4. Discussion: conditions of validity 
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4.1. Scales of the horizontal motions 

One important assumption in our model was that the velocity and time 
scales of horizontal motions due to the shear instability are given by 
equations (5) and (6). This is indeed the general experience with unsta- 
ble shear flows, its physical background being an equilibrium between 
growth and turbulent dissipation of the dominant modes. The growth 
timescale is 

r g ~ l h /Av (27) 
while the dissipation timescale is 

r d ~ ll/v h (28) 

Equating these timescales and substituting ~ l^Vh for the turbulent 
diffusivity, we indeed find Vh ~ Aw. 

But this line of argument is naive for several reasons. First, in the 
strongly stratified, quasi-2D flow in the tachocline the main mecha- 
nism of dissipation is not the nonlinear transfer between 2D modes but 
the viscous dissipation (due to overturning small-scale 3D turbulence) 
between neighbouring slabs. Thus, the dissipation timescale is rather 

r d ~ l\lv. (29) 

Nevertheless, substituting here the formula (7), derived assuming the 
validity of equations (5) and (6), the result is again Vh ~ Au, confirming 
that our assumption was consistent. 

Secondly, equation (27) may not be a valid expression of the growth 
timescale if the instability is only slightly supercritical. Linear 2D hy- 
drodynamical stability analysis indeed shows that the horizontal shear 
in the tachocline is close to marginal stability (Dziembowski and Koso- 
vichev, 1987, Charbonneau, Dikpati, and Gilman, 1999). The weakly 
nonlinear extension of the analysis (Garaud, 2001b) indicates that the 
evolution of unstable modes is such that they modify the mean shear 
profile towards increasing stability, thus ensuring marginal stability, 
low growth rates and low velocity amplitudes. As, however, 3D effects 
and magnetic fields are known to lead to strong instability with much 
weaker shear (Gilman and Dikpati, 2000, Dikpati and Gilman, 2001), 
one may expect that this near-criticality is a feature limited to 2D HD 
models only. 

Somewhat surprisingly, the 3D thin shell simulations of Miesch (2002) 
also show rather low nonaxisymmetric horizontal flow amplitudes, at 
least one order of magnitude below the shear amplitude. This is so 
despite the strongly supercritical rotation profile chosen, which is not 
significantly altered by the feedback of the nonaxisymmetric flow. One 
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may speculate that this is the result of the rather high numerical 
viscosities applied. 

Further study of the flows driven by the horizontal shear instability 
is clearly important. In addition to clarifying the issue of characteristic 
scales of this flow, such studies may shed more light on the question 
of the importance and sense of any direct angular momentum trans- 
port by these motions. At any rate, a lower velocity amplitude or a 
longer timescale for the horizontal flows would result in lower critical 
field strengths for the secondary shear instability, resulting in an even 
thinner tachocline confined by even weaker fields in the framework of 
the present model. 

4.2. Other sources of turbulence 

Another important constraint in the present model was the assumption 
that the secondary shear instability is the only source of turbulence 
in the tachocline. While a number of other instabilities in rotating 
radiative stellar interiors are known (cf. Tassoul, 1978, Zahn, 1983, 
Spruit, 1999), they are mostly acting throughout the radiative interior 
where constraints such as chemical mixing in stellar models impose 
rather strict limits on turbulence. Thus, in the tachocline we need 
"extra" sources of turbulence, for which shear instabilities are the most 
plausible candidate. 

More important than rotational instabilities is the possibility of 
nonadiabatic overshooting convection. The recent solution of the Rey- 
nolds moment equations for the overshoot layer by Marik and Petrovay 
(2003) indicates that the total extent of nonadiabatic overshoot is negli- 
gible. On the other hand, Xiong and Deng (2001) find a very significant 
nonadiabatic overshoot. As both models have their own shortcomings, 
thgis issue is not quite settled yet. At any rate, if there is a significant 
nonadiabatic overshoot, the secondary shear instability would only set 
in below the overshoot layer, where the turbulent diffusivity is reduced 
below the critical value; in this layer, our model could still be considered 
valid. Thus, the main effect of any nonadiabatic overshoot would be 
to "shift" the resulting tachocline to somewhat lower depths (unless 
a stronger dynamo field reduces the horizontal shear already in the 
overshoot region). 



5. Conclusion 

This paper presented the first consistent model for the turbulent tacho- 
cline, with the turbulent diffusivity computed within the model instead 
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of being specified arbitrarily. For the origin of the 3D turbulence a new 
mechanism was proposed, the secondary vertical shear instability of 2D 
motions taking place in thin horizontal slabs. 

Independent evidence for the mechanism proposed here is hard to 
find, as the conditions in the tachocline are extreme by terrestrial stan- 
dards, and direct numerical simulations of stratified shear flows are 
currently limited to much lower values of the Richardson and Reynolds 
numbers. Nevertheless, the presence of a quasi-2D, layered structure, 
and a strongly anisotropical turbulent diffusivity tensor have been con- 
firmed in numerical simulations by Metais and Herring (1989), Kalten- 
bach, Gerz, and Schumann (1994) and Jacobitz and Sarkar (1998). 
Indications for the presence of a secondary shear instability were also 
found in the recent simulations of Werne and Fritts (1999). 

A formula for the diffusivity due to the turbulence generated by the 
secondary shear, equation (9) was derived and applied in a simple ID 
model of the tachocline, taking into account the Maxwell stresses due to 
the dynamo field. As our analysis does not consider the effects of spher- 
ical geometry, rotation, non-adiabatic convective overshoot, meridional 
circulation, or the effect of magnetic fields on stability, it should only 
be regarded as a first step towards a more comprehensive theoretical 
analysis of the problem of turbulence in the solar tachocline. 

A thin tachocline in our models can only be produced if the os- 
cillatory dynamo magnetic field exceeds a critical limit. At this limit, 
the thickness of the tachocline is ~ 5 Mm and the mean value of the 
field strength in the tachocline is 200 G. The total magnetic flux in the 
tachocline in this case agrees well with the total poloidal flux crossing 
the equatorial plane in the poloidal flux transport model of Petrovay 
and Szakaly (1999). 

On the other hand, the tachocline resulting from our model is un- 
comfortably thin when compared to helioseismic constraints. Existing 
seismic calibrations of tachocline thickness (Basu and Antia, 2001) in- 
dicate a significantly thicker shear layer. These calibrations are based 
on the use of simple predefined fitting profiles, characterized by a single 
length scale. In contrast, our model indicates that for slightly subcritical 
magnetic field strengths the profile is more complicated, the fast linear 
cutoff in v within the skin layer being replaced by a shallow but deep 
profile of type (26) below. It remains to be seen, whether such a two- 
tiered -u-profile can yield an equal or better fit to helioseismic data than 
the more conventional profiles. 

The thin tachocline resulting from the present model is also hard 
to reconcile with the gradual depletion of lithium in the atmospheres 
of Sun-like stars during their lifetimes. Lithium is destroyed by nuclear 
reactions in layers below z ~ 40 Mm only, so a mixing characterized 
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by a diffusivity of at least 10 3 cm 2 /s must be present as far down as 
that depth. While our prescription v m = 10 7 cm 2 /s does not allow 
a firm statement on this issue, the very sharp cutoff of the v-curve 
in the right-hand column of Fig. 2 does not seem to indicate that any 
significant level of turbulence could be maintained at such great depths. 
The two-tiered character of our profiles might admit some "leakage" of 
turbulence into much deeper layers, but this would imply some rather 
implausible fine-tuning of the dynamo field strength. 

One obvious shortcoming of the present models is their simplified 
treatment of the time development and of the geometry. The develop- 
ment of more realistic, axially symmetric spherical models employing 
the viscosity formula (9) is in progress. 
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